Differential expression analysis using DESeq2 for matched plasma & serum samples

Match plasma & serum samples

import pandas as pd
samples = pd.read_csv('data/sample_sheet.csv')
# must be Healthy Control Study and must pass MiSeq QC
samples = samples.loc[(samples['Study'] == 'Healthy Controls') & (samples['MISEQ.QC.PASS'] == 'PASS')]  
samples = samples.set_index('MT.Unique.ID').sort_values(by=['Participant.ID', 'Source'])
# capitalize & strip whitespace for consistency
for column in ['Gender', 'Race', 'Source']:
    samples[column] = samples[column].str.capitalize()
    samples[column] = samples[column].str.strip()
# use correct ontology terms
race_ontology = {'Asian': 'Asian',
 'Black or african american': 'African_American',
 'Mixed/asian & white': 'Multiracial',
 'Mixed/asian &black': 'Multiracial',
 'Mixed/black, white, asian': 'Multiracial',
 'Native hawiian or other pacific islander': 'Pacific_Islander',
 'Pacific islander': 'Pacific_Islander',
 'White': 'White'}
for id in samples.index:
    race = samples.at[id, 'Race']
    samples.at[id, 'Race'] = race_ontology[race] if race in race_ontology else 'Multiracial'
# get matched plasma & serum samples
mir_counts = pd.read_csv("data/get_canonical/canon_mir_counts.csv", index_col=0)
samples.index = pd.Index(['X' + str(row) for row in samples.index])
serum_part_ids = set(samples.loc[samples['Source'] =='Serum']['Participant.ID'])
matched_samples = samples.loc[samples['Participant.ID'].isin(serum_part_ids)]
matched_mir_counts = mir_counts[matched_samples.index]
matched_samples.to_csv('data/matched_plasma-serum_samples.csv')
matched_mir_counts.to_csv('data/matched_plasma-serum_mir_counts.csv')

Load the sample data and miRNA counts

samples <- read.csv('data/matched_plasma-serum_samples.csv')
samples <- subset(samples, Library.Generation.Set != "SetRecheck" )  # exclude SetRecheck samples
# sample X11 is outlier -- remove it and its match
outlier_id <- samples[samples$X=="X11",]$Participant.ID
samples <- samples[samples$Participant.ID != outlier_id,]
counts <- read.csv('data/matched_plasma-serum_mir_counts.csv')
samples$Participant.ID <- factor(samples$Participant.ID)  # ID is categorical, not numerical
samples$Sample.ID <- factor(samples$Sample.ID)
rownames(counts) = counts$X
rownames(samples) = samples$X
counts$X <- NULL  # remove extra column
counts <- counts[,rownames(samples)]
head(samples)

head(counts)

Filter

library(edgeR)
design <- model.matrix(~0+Source+Race+Gender+Age, samples)
dge = DGEList(counts = counts, samples = samples)
# require miRNAs to have CPM > 1 in at least 2 samples
countsPerMillion <- edgeR::cpm(dge)
countCheck <- countsPerMillion > 1
head(countCheck)

                 X1  X107    X2  X108   X14  X110   X15   X18   X21  X113   X22  X114   X26  X115   X40  X116   X46  X119   X48  X120   X52
hsa-let-7a-3p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
hsa-let-7a-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
hsa-let-7b-3p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
hsa-let-7b-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
hsa-let-7c-3p FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
hsa-let-7c-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
               X121   X56  X122   X59  X123   X60  X124  X62  X125   X64  X126   X66  X127   X67  X128   X68  X129   X69  X130  X131   X77
hsa-let-7a-3p  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE
hsa-let-7a-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
hsa-let-7b-3p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
hsa-let-7b-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
hsa-let-7c-3p FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
hsa-let-7c-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
               X132  X133   X88  X134   X89  X135   X90  X136  X93  X137  X138
hsa-let-7a-3p FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE TRUE  TRUE  TRUE
hsa-let-7a-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE
hsa-let-7b-3p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE
hsa-let-7b-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE
hsa-let-7c-3p FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
hsa-let-7c-5p  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE TRUE  TRUE  TRUE

keep <- which(rowSums(countCheck) >= 2) 
dge <- dge[keep,]

Explore variance

library(SingleCellExperiment)
library(scater)
reads_sce <- SingleCellExperiment(assays=list(counts=dge$counts),  colData=dge$samples)
# remove unexpressed miRNAs
keep_feature <- rowSums(counts(reads_sce) > 0) > 0
reads_sce <- reads_sce[keep_feature, ]
reads_sce <- calculateQCMetrics(reads_sce)

Note that the names of some metrics have changed, see 'Renamed metrics' in ?calculateQCMetrics.
Old names are currently maintained for back-compatibility, but may be removed in future releases.

# create separate plasma & serum SCEs
plasma_samples <- subset(colData(reads_sce), Source == "Plasma")
plasma_counts <- as.data.frame(counts(reads_sce))[rownames(plasma_samples)]
plasma_sce <- SingleCellExperiment(assays=list(counts=as.matrix(plasma_counts)), colData=plasma_samples)
plasma_sce <- calculateQCMetrics(plasma_sce)

Note that the names of some metrics have changed, see 'Renamed metrics' in ?calculateQCMetrics.
Old names are currently maintained for back-compatibility, but may be removed in future releases.

serum_samples <- subset(colData(reads_sce), Source == "Serum")
serum_counts <- as.data.frame(counts(reads_sce))[rownames(serum_samples)]
serum_sce <- SingleCellExperiment(assays=list(counts=as.matrix(serum_counts)), colData=serum_samples)
serum_sce <- calculateQCMetrics(serum_sce)

Note that the names of some metrics have changed, see 'Renamed metrics' in ?calculateQCMetrics.
Old names are currently maintained for back-compatibility, but may be removed in future releases.

# log transform
cpm(reads_sce) <- calculateCPM(reads_sce)
reads_sce <- normalize(reads_sce)

using library sizes as size factors

logcounts(reads_sce) <- log2(calculateCPM(reads_sce) + 1)
# visualize
hist(reads_sce$total_counts, breaks=100)  # counts per sample

hist(reads_sce$total_features, breaks=100)  # counts per miRNA

plotQC(reads_sce, type = "highest-expression") + ggtitle("Plasma & Serum")

plotQC(plasma_sce, type = "highest-expression") + ggtitle("Plasma")

plotQC(serum_sce, type = "highest-expression") + ggtitle("Serum")

plotQC(reads_sce, type="explanatory-variables", variables=c("Index", "Participant.ID", "Collection.Date", "Library.Generation.Set", "MiSeq.QC.Run", "Source", 'total_features', "Age", "Race", "Sex"))

variable Sex not found in colData(object).
                     Please make sure colData(object)[, variable] exists. This variable will not be plotted.

Examine sources of variation

plotPCA(reads_sce, exprs_values = "logcounts", colour_by = "Library.Generation.Set", size_by = "total_features", shape_by = "Source")

non-plotting arguments like 'exprs_values' should go in 'run_args'

for (var in c("total_features", "Age", "Library.Generation.Set", "Index", "Participant.ID", "Source", "Race", "Gender", "Collection.Date")) {
  print(
    plotQC(reads_sce, type = "find-pcs", exprs_values = "logcounts", variable = var)
    )  + ggtitle(var) 
  }

Normalize & Remove unwanted sources of variation

library(RUVSeq)
library(ggplot2)
library(mvoutlier)
dge <- calcNormFactors(dge)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
lrt <- glmLRT(fit, coef=2)
res <- residuals(fit, type="deviance")
set <- newSeqExpressionSet(dge$counts, phenoData=dge$samples)
for(k in 1:10) {
  ruvr_set <- RUVr(set, row.names(dge), k=k, res)
  assay(reads_sce, paste("RUVr k=", toString(k))) <- log2(t(t(assayData(ruvr_set)$normalizedCounts) / colSums(assayData(ruvr_set)$normalizedCounts) * 1e6) + 1)
}
for(n in assayNames(reads_sce)) {
  print(
        plotPCA(
            reads_sce,
            colour_by = "Library.Generation.Set",
            size_by = "total_features",
            shape_by = "Source",
            exprs_values = n
        ) + ggtitle(n)
  )
}

non-plotting arguments like 'exprs_values' should go in 'run_args'

Detect outliers

reads_sce <- runPCA(reads_sce, use_coldata = TRUE, detect_outliers = TRUE)

failed to find 'pct_counts_feature_control' in column metadatafailed to find 'total_features_by_counts_feature_control' in column metadatafailed to find 'log10_total_counts_endogenous' in column metadatafailed to find 'log10_total_counts_feature_control' in column metadata

outliers <- colnames(reads_sce)[reads_sce$outlier]
head(outliers)

character(0)

Examine sources of variance after removing unwanted variation

for (var in c("total_features", "Age", "Library.Generation.Set", "Index", "Participant.ID", "Source", "Race", "Gender", "Collection.Date")) {
  print(
    plotQC(reads_sce, type = "find-pcs", exprs_values = "RUVr k= 2", variable = var)
    )  + ggtitle(var) 
}

Visualize top plasma-vs-serum miRNAs

Visualize EdgeR results

topTags(lrt)
# TODO: filter by p-value and/or FDR
# TODO: heatmap

Run the DESeq2 pipeline & Visualize results

library(DESeq2)  # TODO: use filtered & RUVSeq-corrected data
deseq2Data <- DESeqDataSetFromMatrix(countData = counts, design = ~ Participant.ID + Source, colData = samples)
deseq2Data <- DESeq(deseq2Data)
results <- results(deseq2Data, cooksCutoff=FALSE, independentFiltering=FALSE)
head(as.data.frame(results))

Create a heatmap of the variance-stabilizing-transformed expression data

library(pheatmap)   # TODO: only use top DE miRNAs
vstData <- varianceStabilizingTransformation(deseq2Data, blind=FALSE)  # transform the data
select <- order(rowMeans(counts(deseq2Data,normalized=TRUE)), decreasing=TRUE)
Source <- colData(deseq2Data)[,c("Source")]
annotations <- as.data.frame(Source)
rownames(annotations) <- colnames(deseq2Data)
pheatmap(assay(vstData)[select,], cluster_rows=TRUE, show_rownames=FALSE, cluster_cols=TRUE, annotation_col=annotations)

Create a heatmap of sample-to-sample distances

sampleDists <- dist(t(assay(vstData)))  # compute sample-to-sample distances
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vstData$Participant.ID, vstData$Source, sep=" - ")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)

Create a PCA plot

pcaData <- plotPCA(vstData, intgroup=c("Source"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=Source)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()

save.image("diff_expr_plasma_vs_serum.RData")